function dS = NBody(~,S)
    
    u = [132712440041.93938; 22032.09; 324858.63; 398600.44; 42828.3; 0];
    %[km^3/s^2]Gravitational parameters of the Sun, Mercury, Venus, the Earth, and Mars.
    
    N = numel(u);
    %[]Number of bodies that will be analyzed.
    
    R = reshape(S(1:3 * N),3,N);
    %[m]Position matrix.
    
    A = zeros(3,N);
    %[m/s^2]Allocates memory for the acceleration matrix.
    
    for i = 1:N
        %[]The number of columns in the position matrix is the number of bodies being analyzed.
        
        for k = 1:N
            %[]For each body being analyzed, the remaining bodies will perturb its motion.
            
            if i == k
                
                continue;
                
            end
            %[]The body in question does not perturb itself.
            
            A(:,i) = A(:,i) + u(k) * (R(:,k) - R(:,i)) / norm(R(:,k) - R(:,i))^3;
            %[m/s^2]Updates the acceleration vector of the current body.
            
        end
        
    end
    
    dS = zeros(6 * N,1);
    %[m/s,m/s^2]Allocates memory for the differential state vector.
    
    dS(1:3 * N) = S(3 * N + 1:6 * N);
    %[m/s]Velocity components.
    
    dS(3 * N + 1:6 * N) = reshape(A,3 * N,1);
    %[m/s^2]Acceleration components.
    
end
%===================================================================================================